#include "vector.h"
#include "matriz.h"
#include "sistemas.h"

/* Probar las funciones de sistemas.cpp */
/* que son:

void invertirTriangInf(Matriz& A, Matriz& res, int n); tick
void matrizPorVector(Matriz &A, Vector &x, Vector& res, int n); tick
void matrizPorDiagonal(Matriz& A, Vector &d, Vector& res, int n); tick
void diagonalPorMatriz(Vector& d, Matriz& A, Vector& res, int n); tick
void backwardSubs(Matriz& L, Vector& b, Vector& x, int n); tick
void forwardSubs(Matriz& U, Vector& b, Vector& x, int n); tick
void PLUPivoteo(const Matriz& M, Vector& b, Vector& x, int n);

*/

int main()
{
	cout.precision(4);
	int n = 5;
	
	Matriz A (n, Vector(n, 0)); /* L */
	Matriz B (n, Vector(n, 0)); /* U */
	
	forn(i, n)
	{
		//~ A[i][i] = 1;
		
		fornd(j, i+1)
		{
			A[i][j] = (i + j + 1);
			//~ A[i][j] = (i + j + 1)/5.42;
		}
	}
	
	forn(i, n)
	{
		fornx(j, i, n)
		{			
			B[i][j] = i*i + j + 1;
		}
	}
		
	cout << "A" << endl;
	mostrarMatrizAlpha(A, n);
	cout << "B" << endl;	
	mostrarMatriz(B, n);
	Matriz invA (n, Vector(n));
	Matriz invB (n, Vector(n));
	invertirTriangInf(A, invA, n);	
	cout << "Inversa de A" << endl;	
	mostrarMatriz(invA, n);
	//~ cout << "Inversa de B" << endl;	
	//~ mostrarMatriz(invB, n);	
	
	Matriz resA (n, Vector(n));
	mul(A, invA, resA, n);
	cout << "A * inversa de A" << endl;
	mostrarMatriz(resA, n);
	
	n = 4;
	Matriz C (n, Vector(n, 0));
	Vector x (n, 0);
	Vector res (n, 0);
	
	C[0][0] = 0; C[0][1] = 1; C[0][2] = 5; C[0][3] = 7;
	C[1][0] = 15; C[1][1] = 2; C[1][2] = 3; C[1][3] = 8;
	C[2][0] = 6; C[2][1] = 4; C[2][2] = 1; C[2][3] = 10;
	C[3][0] = 2; C[3][1] = 5; C[3][2] = 2; C[3][3] = 7;
	
	x[0] = 4; x[1] = 9; x[2] = 11; x[3] = 1;
	
	matrizPorVector(C, x, res, n);
	cout << "C" << endl;
	mostrarMatrizAlpha(C, n);
	cout << "x" << endl;
	mostrarVectorAlpha(x, n);
	cout << "C*x" << endl;
	mostrarVectorAlpha(res, n);
	
	Matriz resMatriz (n, Vector(n));
	matrizPorDiagonal(C, x, resMatriz, n);
	cout << "C*x (con x la diagonal de una matríz)" << endl;
	mostrarMatrizAlpha(resMatriz, n);
	diagonalPorMatriz(x, C, resMatriz, n);
	cout << "x*C (con x la diagonal de una matríz)" << endl;
	mostrarMatrizAlpha(resMatriz, n);
	
	/* Ahora en double. */	
	forn(i, n)
		forn(j, n)
			C[i][j] /= 5.42;
			
	forn(i, n)
		x[i] /= 5.42;
		
	matrizPorVector(C, x, res, n);
	cout << "C" << endl;
	mostrarMatrizAlpha(C, n);
	cout << "x" << endl;
	mostrarVectorAlpha(x, n);
	cout << "C*x" << endl;
	mostrarVectorAlpha(res, n);
	matrizPorDiagonal(C, x, resMatriz, n);
	cout << "C*x (con x la diagonal de una matríz)" << endl;
	mostrarMatrizAlpha(resMatriz, n);
	diagonalPorMatriz(x, C, resMatriz, n);
	cout << "x*C (con x la diagonal de una matríz)" << endl;
	mostrarMatrizAlpha(resMatriz, n);
	
	n = 5;
	
	//~ forn(i, n-1)
		//~ x[i] *= 5.42;
	
	x.push_back(2);
	res.push_back(0);
	matrizPorVector(A, x, res, n);
	
	cout << "L" << endl;
	mostrarMatrizAlpha(A, n);
	cout << "x" << endl;
	mostrarVectorAlpha(x, n);
	cout << "L*x" << endl;
	mostrarVectorAlpha(res, n);
		
	backwardSubs(A, res, x, n);	
	
	cout << "x?" << endl;
	mostrarVectorAlpha(x, n);	
	cout << endl;
	
	matrizPorVector(B, x, res, n);
	cout << "U" << endl;
	mostrarMatrizAlpha(B, n);
	cout << "x" << endl;
	mostrarVectorAlpha(x, n);
	cout << "U*x" << endl;
	mostrarVectorAlpha(res, n);
	
	forwardSubs(B, res, x, n);
	
	cout << "x?" << endl;
	mostrarVectorAlpha(x, n);
	
	n = 4;

	Matriz L (n, Vector(n, 0));
	Matriz U (n, Vector(n, 0));
	
	forn(i, n)
	{
		L[i][i] = 1;
		U[i][i] = i + 1;
	}
		
	L[1][0] = 2; L[2][0] = 3; L[2][1] = 2;
	L[3][0] = 4; L[3][1] = 3; L[3][2] = 2;
	
	U[0][1] = 5; U[0][2] = 8; U[0][3] = 10;
	U[1][2] = 7; U[1][3] = 6; U[2][3] = 9;
	
	cout << "L" << endl;
	mostrarMatrizAlpha(L, n);
	cout << "U" << endl;
	mostrarMatrizAlpha(U, n);
	
	Matriz LUM (n, Vector(n));	
	mul(L, U, LUM, n);
	
	//~ x[0] = 1; x[1] = 1; x[2] = 1; x[3] = 1;
	//~ forn(i, n)
		//~ x[i] *= 5.42;
	
	matrizPorVector(LUM, x, res, n);	
	
	cout << "LU" << endl;
	mostrarMatrizAlpha(LUM, n);
	cout << "x" << endl;
	mostrarVectorAlpha(x, n);
	cout << "LU*x" << endl;
	mostrarVectorAlpha(res, n);
	
	Matriz resLU (n, Vector(n));
	
	PLUPivoteo(LUM, res, x, n);
	//~ PLUPivoteo(LUM, res, x, resLU, n);
	//~ cout << "res LU" << endl;
	//~ mostrarMatrizAlpha(resLU, n);
	//~ 
	//~ Matriz LPrima (n, Vector(n, 0));
	//~ Matriz UPrima (n, Vector(n, 0));
	//~ 
	//~ forn(i, n)
	//~ {
		//~ forn(j, i)
		//~ {
			//~ LPrima[i][j] = resLU[i][j];
		//~ }
		//~ 
		//~ LPrima[i][i] = 1;
		//~ 
		//~ fornx(j, i, n)
		//~ {
			//~ UPrima[i][j] = resLU[i][j];
		//~ }
	//~ }
	//~ 
	//~ cout << "LPrima" << endl;
	//~ mostrarMatrizAlpha(LPrima, n);
	//~ cout << "UPrima" << endl;
	//~ mostrarMatrizAlpha(UPrima, n);
	//~ 
	//~ mul(LPrima, UPrima, resLU, n);
	//~ 
	//~ cout << "LPrima * UPrima" << endl;
	//~ mostrarMatrizAlpha(resLU, n);
	cout << "x?" << endl;
	mostrarVector(x, n);
		
	//~ n = 3;
	//~ Matriz pruebaLU (n, Vector(n));
	//~ 
	//~ pruebaLU[0][0] = 8.0; pruebaLU[0][1] = 4.0; pruebaLU[0][2] = 5.0;
	//~ pruebaLU[1][0] = 2.0; pruebaLU[1][1] = 3.0; pruebaLU[1][2] = 7.0;
	//~ pruebaLU[2][0] = 1.0; pruebaLU[2][1] = 2.0; pruebaLU[2][2] = 4.0;
	//~ LU(pruebaLU, res, x, resLU, n);
	
	//~ n = 4;
	//~ Matriz Lower (n, Vector(n, 0));
	//~ Matriz inversaLower (n, Vector(n));
	//~ 
	//~ forn(i, n)
		//~ Lower[i][i] = 1;
		//~ 
	//~ Lower[1][0] = 5; Lower[2][0] = 27; Lower[2][1] = 9;
	//~ Lower[3][0] = 31; Lower[3][1] = 11; Lower[3][2] = 2;
	//~ 
	//~ invertirTriangInf(Lower, inversaLower, n);
	//~ 
	//~ mostrarMatriz(inversaLower, n);
	//~ 
	//~ Matriz identidad (n, Vector(n));	
	//~ mul(Lower, inversaLower, identidad, n);	
	//~ mostrarMatriz(identidad, n);
	
	n = 4;
	Matriz M (n, Vector(n));
	
	M[0][0] = 5; M[0][1] = 27; M[0][2] = 8; M[0][3] = 21;
	M[1][0] = 9; M[1][1] = 15; M[1][2] = 10; M[1][3] = 34;
	M[2][0] = 2; M[2][1] = 14; M[2][2] = 13; M[2][3] = 42;
	M[3][0] = 17; M[3][1] = 3; M[3][2] = 2; M[3][3] = 1;
	
	x[0] = 24; x[1] = 15; x[2] = 197; x[3] = 205;
	
	/* Ahora en double. */	
	forn(i, n)
		forn(j, n)
			M[i][j] /= 5.42;
			
	forn(i, n)
		x[i] /= 5.42;
		
	cout << "M" << endl;
	mostrarMatrizAlpha(M, n);
	cout << "x: ";
	mostrarVectorAlpha(x, n);
	
	matrizPorVector(M, x, res, n);
	
	cout << "res: ";
	mostrarVectorAlpha(res, n);
	
	PLUPivoteo(M, res, x, n);
	
	cout << "x?: ";
	mostrarVectorAlpha(x, n);
	
	Matriz Hil (n, Vector(n));
	Hilbert(Hil, n);
	
	cout << "res: ";
	mostrarVector(res, n);
	
	buscarMatriz(M, Hil, res, x, n);
	cout << "Matriz" << endl;
	mostrarMatriz(M, n);
	
	matrizPorVector(M, x, res, n);
	
	cout << "res? ";
	mostrarVector(res, n);
	
	return 0;
}
